function [df_sim, dfull_sim, coefs_sim] = data_simulate( ...
    xi,       ...
    eta,      ...
    coefs,    ...
    df,       ...
    dfull,    ...
    nmarkets, ...
    fe_miss,  ...
    version,  ...
    seed      ...
)
    
    markets  = unique(df.cdid);
    nmarket  = nmarkets;

    xz_market = datasample(markets,nmarket);
    ew_market = datasample(markets,nmarket);

    % store derivative matrix
    der_sim      = zeros(max(df.prodid),max(df.prodid),nmarket);
    der_mod      = zeros(max(df.prodid),max(df.prodid),nmarket);

    % creagte shelf variable
    brndcd    = findgroups(df.brndid,df.cdid);
    shelf     = accumarray(brndcd,1);
    df.shelf  = shelf(brndcd);

    
    % UNIFORM INCOME DISTRIBUTION ACROSS MARKETS AND NORMALIZE AVG INCOME
    C                 = size(dfull,2); 
    cityyearid    = grp2idx(df.yearid*100+df.cityid);
    mean_inc      = accumarray(cityyearid,df.minc,[],@mean);
    income_draws  = cell2mat(arrayfun(@(x) accumarray(cityyearid,dfull(:,x),[],@(x)mean(x)),1:C,'un',0))+mean(mean_inc);
    income_norm   = income_draws./mean_inc;
    income_norm   = income_norm(:);
    income_sample = datasample(income_norm,C)';
    minc_new      = coarsen_income(mean_inc);
    df.minc       = minc_new(cityyearid);

    
    % DRAW (X^SIM,Z) WITH REPLACEMENT

    df_sim = table();
    for m = 1:nmarket
        select    = (df.cdid==xz_market(m));
        cdid_sim  = array2table(m*ones(sum(select),1),'VariableNames',{'cdid_sim'});
        df_sim    = [df_sim;df(select, :),cdid_sim];
    end
    df_sim.marketsize = round(df_sim.q_jt./df_sim.s_jt);
    dfull_sim         = repmat(income_sample,size(df_sim.id2)).*df_sim.minc;
    dfull_sim         = dfull_sim - mean(mean(dfull_sim));


    % DRAW (XI,OMEGA) WITH REPLACEMENT

    J = max(df.prodid);
    ew_final = [];
    for m = 1:nmarket
        select = (df.cdid==ew_market(m));
        ew_sim = [xi(select),eta(select)];
        prod_sim = df.prodid(select);
        ew_new = nan(J,2);
        ew_new(prod_sim,:) = ew_sim;
        ew_fill = datasample([xi,eta],J);
        ew_new(isnan(ew_new(:,1)),:) = ew_fill(isnan(ew_new(:,1)),:);
        ew_final = [ew_final;ew_new(df_sim.prodid(df_sim.cdid_sim==m),:)];
    end
    df_sim.xi = ew_final(:,1);
    df_sim.eta = ew_final(:,2);


    % COLLAPSE DATE FIXED EFFECTS INTO MONTH FIXED EFFECTS

    sprod    = f_selectfe(df.prodid, df.prodid, 0);
    scity    = f_selectfe(df.cityid, df.cityid, 1);
    sdate    = f_selectfe(df.dateid, df.dateid, 1);
    sdatem   = accumarray(df.dateid,df.montid,[],@mean); 
    montD_fe = accumarray(sdatem(2:end),coefs.sigma_D(sdate+max(df.prodid)-1),[],@mean);
    montS_fe = accumarray(sdatem(2:end),coefs.sigma_S(sdate+max(df.cityid)+max(df.prodid)-2),[],@mean);


    % COLLAPSE PRODUCT FIXED EFFECTS INTO BRAND FIXED EFFECTS

    sbrand   = accumarray(df.prodid,df.brndid,[],@mean); 
    brndD_fe = accumarray(sbrand,coefs.sigma_D(sprod),[],@mean);
    brndD_fe = brndD_fe(sbrand);

    % CALCULATE FIXED EFFECTS, MARGINAL COST, ETC

    % Fixed effects
    montfe = df_sim.montid == unique(df_sim.montid)';
    datefe = df_sim.dateid == unique(df_sim.dateid)';
    prodfe = df_sim.prodid == unique(df_sim.prodid)';
    cityfe = df_sim.cityid == unique(df_sim.cityid)';
    fesd   = [prodfe montfe];
    fess   = [prodfe cityfe(:, 2:end) montfe];

    % Variables for demand and supply
    x   = [df_sim.p_jt fesd];
    if strcmp(version,'_sh')
        x = [x,df_sim.shelf];
    end
    w   = [df_sim.mcpost df_sim.dist fess];
    xnonlin = [df_sim.p_jt df_sim.cons df_sim.calor];

    % other parameters
    
    prodD_fe = (1-fe_miss)*coefs.sigma_D([sprod]) + fe_miss*brndD_fe;
    sigma_S  = [coefs.sigma_S([sprod; scity+max(df.prodid)-1]);montS_fe];
    sigma_D  = [prodD_fe;montD_fe];
    if strcmp(version,'_sh')
        sigma_D = [sigma_D;coefs.sigma_D(end)];
    end
    theta    = [coefs.alpha; sigma_D];
    bc       = (df_sim.yearid > 4) .* coefs.kappa;
    mcost    = w * [coefs.gamma; sigma_S] + df_sim.eta;

    % make sure all costs are ge 0
    mcost(mcost<0) = 0;

    % FIND EQUILIBRIUM 

    df_sim.s_jt = zeros(size(mcost));
    df_sim.p_jt = zeros(size(mcost));

    tic;
    options = optimoptions('fsolve', 'Display', 'none');
    for m = 1:nmarket
        if mod(m, 10) == 0
            fprintf('Time passed: %.2f minutes', toc/60);
            disp(' ')
            fprintf('Progress: %4d / %4d', m, nmarket);
            disp(' ')
            disp(' ')
        end
        sel = (df_sim.cdid_sim == m);
        foc = @(p) f_equilibrium(  ...
            p,                     ...
            1,                     ...
            theta,                 ...
            coefs.Pi,              ...
            coefs.rho,             ...
            x(sel, :),             ...
            xnonlin(sel, :),       ...
            df_sim.xi(sel, :),     ...
            dfull_sim(sel, :),     ...
            mcost(sel, :),         ...
            unique(bc(sel)),       ...
            df_sim.bmc(sel),       ...
            df_sim.firmid(sel)     ...
        );

        price_jt = fsolve(foc, df_sim.p_jt(sel), options);
        [~, share_jt, ~, dshare_jt, ~, der_jt] = foc(price_jt);

        df_sim.p_jt(sel)    = price_jt;
        df_sim.s_jt(sel)    = share_jt;
        df_sim.dsdp_jt(sel) = dshare_jt;

        der_model                  =  share_jt*share_jt';
        der_model(logical(eye(length(share_jt)))) = -share_jt.*(1-share_jt);

        der_sim(df_sim.prodid(sel),df_sim.prodid(sel),m)      = der_jt;
        der_mod(df_sim.prodid(sel),df_sim.prodid(sel),m)      = der_model;
    end

    % PREPARE OUTPUT

    coefs_sim          = coefs;
    coefs_sim.sigma_D  = sigma_D;
    coefs_sim.sigma_S  = sigma_S;
    coefs_sim.dparams  = [coefs.alpha; coefs.Pi; coefs.rho];
    coefs_sim.sparams  = [coefs.kappa; coefs.gamma];
    coefs_sim.der_sim  = der_sim;
    if seed == 1
        coefs_sim.J       = J;
        coefs_sim.der_mod = der_mod;

    df_sim = removevars(df_sim,{'q_jt',    ...
                                'inshr',   ...
                                'outshr',  ...
                                's_jgt',   ...
                                'logodds', ...
                                'cons',    ...
                                });

end
